Numerical calculation of the combinatorial entropy of partially ordered ice 
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Using a one-parameter case as an example, we demonstrate that multicanonical simulations allow for 
accurate estimates of the residual combinatorial entropy of partially ordered ice. For the considered 
case corrections to an (approximate) analytical formula are found to be small, never exceeding 0.5%. 
The method allows one as well to calculate combinatorial entropies for many other systems. 



I. INTRODUCTION 

After the discovery of the hydrogen bond it was rec- 
ognized that the unusual properties of water and ice owe 
their existence to a combination of strong directional po- 
lar interactions and a network of specifically arranged 
hydrogen bonds^"^. By experimental discovery^ it was 
found that ice I (ordinary ice) has in the zero temperature 
limit^ a residual entropy S — k In(W^i) > where Wi is 
the number of configurations per molecule. Subsequently 
Linus Pauling^ based the estimate VF^^^"''"^ — 3/2 on the 
ice rules: 

1. There is one hydrogen atom on each bond (then 
called hydrogen bond) . 

2. There are two hydrogen atoms near each oxy- 
gen atom (these three atoms constitute a water 
molecule). 

Pauling's combinatorial estimate turned out to be in ex- 
cellent agreement with subsequent refined experimental 
measurements*. This may be a reason, why it took 
25 years until Onsager and Dupuis^ pointed out that 
Wi = 1.5 is only a lower bound, because Pauling's ar- 
guments for disordered ice omits correlations induced by 
closed loops which are encountered when one requires ful- 
fillment of the ice rules for all molecules. Subsequently 
Nagle^" used a series expansion method to derive the es- 
timate W^f ^"8'° = 1.50685 (15), where the error bar is not 
statistical but reflects higher order corrections of the ex- 
pansion, which are not rigorously under control. 

Groundstate entropy calculations by means of mul- 
ticanonical (MUCA)^^ Markov chain Monte Carlo 
(MCMC) simulations were pioneered by Berg and 
Celik^^. In a recent paper^^ it was shown that this ap- 
proach allows rather easily for an accurate finite-size scal- 
ing estimate of the residual entropy of ice I, Wf^^'~^^ — 
1.50738 (16), where the error bar is now purely statisti- 
cal. In view of eventual higher order finite size correc- 
tions, which are not included in the MUCA error bar, 
there is satisfactory agreement with Nagle^". 

With the advent of neutron scattering technology, it 
became possible to measure the actual hydrogen arrange- 
ments. Besides fully ordered and disordered ice phases. 



there is also evidence for partially ordered ice'^ . 
Based on theoretical groundwork laid by Takagi^^ and 
Minagawa^*, an extension of Pauling's results to partially 
ordered ice was derived bu Howe and Whitworth^^ and 
greatly generalized by MacDowell et al.^°. Comparisons 
with neutron scattering results are also made in Ref.^*^. 
Besides, the combinatorial residual entropy needs to be 
taken into account when one considers the phases of sim- 
ple models for water/ice^^. 

As for disordered ice in Pauling's work, correlations are 
neglected in the analytical estimates^*~^° of the residual 
entropy of partially ordered ice. The magnitude of cor- 
rections is largely unknown. For instance, before the pa- 
per by Howe and Whitworth an erroneous equation was 
used, which was off by up to more than 50% for the en- 
tropy per molecule. Nagle's method appears to be too 
complicated for these situations. In this article we gen- 
eralize the MUCA approach of Ref.^^ to include partial 
order and calculate numerical corrections to the formula 
of Howe and Whitworth^^. Our method is presented in 
section II, details of our numerical implementation are 
given in section HI, followed by the entropy estimates in 
section IV. Summary and conclusion with an outlook on 
other applications are given in the final section V. 

II. THE METHOD AND PRELIMINARIES 

As in^^ we confine our interest to the hexagonal crystal 
structure of which the z = layer is is depicted in Fig. 1. 
Each oxygen atom is located at the center of a tetrahe- 
dron and straight lines (bonds) through the sites of the 
tetrahedron point towards four nearest-neighbor oxygen 
atoms. Distances in this figure are given in units of a lat- 
tice constant a (a = 1 in the figure(, which is chosen to 
be the edge length of the tetrahedra. The distance from 
the center of a tetrahedron to one of its sites is y/3/8a 
and, hence, the oxygen-oxygen distance is y^3/2a. 

This is not the conventional crystallographic definition, 
but convenient for setting up the computer program (see 
below). For each molecule shown one of the surface tri- 
angles of its tetrahedron is placed in the a;y-plane. The 
molecules labeled by u (up) are then at z = a/\/24 above. 
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FIG. 1. (Color online) Lattice structure of the z = layer 
of ice I. The up (u) sites are at z — l/v^ and the down (d) 
sites at z = —1/^/24. For each site three of its four bonds 
to nearest neighbor sites are shown. The fourth bond (to the 
next layers) is in up direction for up and in down direction 
for down sites. 

and the molecules labeled by d (down) at z = —a/^/24: 
below the xy-plane, at the centers of their tetrahedra. 

We define an ordered reference configuration, which 
fulfills the ice rules, by arranging the hydrogen atoms on 
the bonds in the following way: 

1. For z — iz'ia/V6 and even (as shown in Fig. 1 
for iz = 0): For the up oxygens put the hydrogens 
on bonds 2 and 4, for the down oxygens put them 
on bonds 1 and 2. 

2. For z = izia/^/d and iz odd (as shown in Fig. 2 for 
iz — 1): For the up oxygens put the hydrogens on 
bonds 3 and 4, for the down oxygens put them on 
bonds 1 and 3. 

This serves as our ordered reference configuration and we 
denote the hydrogen positions in this configuration by rf,. 

FoUowing^^ we denote the fraction of Hydrogen posi- 
tions, which agree with the reference configuration by /. 
The analytical approximation^^ for the residual entropy 
of configurations which fulfill the ice rules is given by 

]5P(/-p)2(/-p)(l+p-2/)(i+P-2/) ^ ' 

with p = / - 2(1 - v/3/2 - 3/+ l)/3. The probabihty 
that the position of a hydrogen atom agrees by chance 
with the one in the reference configuration is 1/2 for 
disordered ice and Eq. (1) reproduces Pauling's result 
for this case, W^{l/2) = 1.5. For / = 1 all hydro- 
gen positions agree with the reference configuration, and 
W{'(1) — 1. Correlations due to closed loops of hydro- 
gen bonds are neglected in the arguments, which lead to 
Eq. (1). Here they are included numerically. 
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FIG. 2. (Color online) Lattice structure of the next, 
2 = 2 + 2/24, layer of ice I, above the one of Fig. 1. 



The residual entropy of ice I was calculated in^^ by 
performing MUCA simulations for two discrete statistical 
model, which were constructed to satisfy the following 
properties [/3 = l/(fcT)]: 

1. Their total number states (as sampled at /9 = 0) is 
known. 

2. Generically each model fulfills one of the ice rules, 
but not the other. 

3. In their energy groundstates (reached at large 
enough /3) each model fulfills both ice rules. 

The model, which fulfills ice rule 2 generically is called 
6-statc H2O molecule model and has for N molecules 
a total number of 6^ states. The model, which fulfills 
ice rule 1 generically is called 2-state H-bond model and 
has 2^^ = 4^ states. Both systems have similarities 
with Potts models, so that the lattice labeling outlined by 
Fig. 1 and 2 allows one^^ to employ simulation methods 
entirely analogue to those outlined for Potts models in 
Ref.^'*. Groundstate entropy estimates with the 2-state 
H-bond model turned out to be more efficient than those 
with the 6-state H2O molecule model, apparently because 
4 is closer to 1.5 than 6. So we confine our generalization 
for partially ordered ice to the 2-state model. To do the 
same for the 6-state model is straightforward, but the 
simulations are expected to be less efficient. 

In the 2-state H-bond model^'^ we allow two positions 
for each hydrogen nucleus on its bond (close toe either 
one of the two oxygen atoms, which are connected by the 
bond). The energy is defined by 

E=-Y.je{,s,h\,hl,hl,h\), (2) 

where the sum is over all sites (oxygen atoms) of the 
lattice and the function /e is given by 
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fAs,bl,blblbt) = 

{2 for two hydrogen nuclei close to s, 
1 for one or three hydrogen nuclei close to s, 
for zero or four hydrogen nuclei close to s. 

We consider now an additional term 



(3) 



(4) 



TABLE I. Simulation statistics overview. 



iV 




Uy 




Statistics 


Additional /lo-valucs used 


128 


4 


8 


4 


32 X 10" 


0.70 


360 


5 


12 


6 


32 X 10" 


0.65, 0.66, 0.67, 070 


576 


6 


12 


8 


32 X 3 10" 


0.65, 0.66, 0.67, 070 


896 


7 


16 


8 


32 X 9 10" 


0.65, 0.66, 0.67, 0.68, 070 


1600 


8 


20 


10 


32 X 32 10" 


0.66, 0.67, 0.68 



which is the overlap of the actual positions Xb of the 
hydrogen atoms on the bonds h with the reference posi- 
tions ffe . The canonical ensemble of the extended model 
is defined by the Gibbs-Boltzmann weights 



exp(-/Ji; + /iQ). 



(5) 



The coupling parameter h plays pretty much the same 
role as an external magnetic field does for the Ising 
model. 

At /? = the expectation value of the overlap per link 
is readily computed to be^^ 



(g)/3=o = (g)o = (Q)o/(2A^) 



1 



(6) 



and the number of states for which the positions of K 
hydrogen atoms agree with those in the reference config- 
uration is given by the binomial factor 



B{2N,K) 



2N 
K 



(27V)! 



{2N -K)\ K\ 



(7) 



The fraction of correct bonds with respect to the refer- 
ence configuration is given by / = K/{2N). For K « 
(Q)o there will be sufficient statistics so that reweighting 
of the simulation to /3 = can be used to normalize the 
spectral density via the binomial distribution (7). For 
that purpose it is convenient to choose h so that {Q)q 
becomes an integer. Assuming that this is done, we take 
K = {Q)o in the following. 

Using a MUCA weight function 



^MUCA ^ ghQ ^MUCA^£;^ 



(8) 



we can connect the (3 = region, for which the numbers 



of states arc known, to the groundstate Eg 



-2N, for 



which both ice rules are satisfied, and estimate the num- 
ber of states n{Q,Eg) for Q values encountered in the 
groundstates with sufficient statistics in H^^'~^^{Q, Eg) 
by reweighting: 



H 



MUCA 



niQ,Eg) 
B{2N,K) J2EH^^'^^iK,E)/'> 



iQ,Eg)/wr^^iQ,Eg) 



,,MUCA 



{K,E) 



(9) 



Here H {Q, Eg) is the overlap histogram sampled by 
the multicanonical updating in the groundstate ensemble 
and ijJ^ucA f^j^^ j^j ^j-^^ energy histogram sampled for 
the fixed value Q = K. The reweighting is to /3 = with 



h tmchangcd. As the MUCA weights (8) factorize, the 
storage requirements are of order N (not N"^). 

To obtain a working estimate (see chapter 5.1 of^**) of 
the MUCA weights we use the Wang-Landau recursion^^ 
as explained in the next section. The numerical quanti- 
ties encountered in Eq. (9) are often so large that they 
are not allowed by a conventional programing language 
like Fortran" 77 in Real*8 precision. This is overcome by 
using consistently logarithmic coding for which technical 
detail are explained in^"'. 

The actually covered Q range in the groundstate en- 
semble depends on h. Increasing h will shift the range to 
higher Q values. Doing so in small steps, and repeating 
the simulation each time, 



Wi{f) = ^Hn{Q,Eg)] 
is obtained for all desired values of / = Q/{2N). 



(10) 



III. NUMERICAL IMPLEMENTATION 

Using periodic boundary conditions (BCs), our sim- 
ulations are based on the lattice construction of Fig. 1 
and 2. Following closely the method otitlined in chap- 
ter 3.1.1 of^* four index pointers from each molecule to 
the array positions of its nearest neighbor molecules are 
constructed along the directions of the bonds as out- 
lined in Fig. 1. The lattice contains then N = UxUyUz 
molecules, where n^, %, and riz are the numbers of sites 
along the x, y, and z axes, respectively; i^; = 0, . . . , Ux — l, 
iy = 0, . . . , rij, — 1, and = 0, . . . , — 1. The pe- 
riodic BCs restrict the allowed values of Hx, Uy, and 
riz to Hx = 1, 2, 3, . . ., Uy = 4, 8, 12, . . ., and Hz = 
2, 4, 6, ... . Otherwise the geometry does not close prop- 
erly. With the inter-site distance rpo = 2.764 A from 
Ref.^, the physical size of the box is obtained by putting 
the lattice constant to a = 2.257 A, and the physical di- 
mensions of the box are calculated to be Bx = 2nx a, 
By = {uy V3/2) a, Bz = (n^ 4/^/6) a. In our choices of 
nx 

tions at symmetrically sized boxes. 

The h values at which the simulations are performed 
are determined from initially proposed ho values in the 
following way: Prom (6) we calculate {Q)o{ho) and de- 
termine the closest integer K. Then the relation (6) 
is inverted to find the value h for which the relation 
K = {Q)o{h) holds. All our simulations use ho = 0.1, 



vx, riy, and values we aim within reasonable limita- 
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0.2, 0.3, 0.4, 0.5, and 0.6. Additional /lo-values are listed 
in Tabic I, which gives an overview of our lattices and 
MUCA production statistics. The statistics is in sweeps 
(i.e., updates per molecule) and repeated 32 times. For 
each of the 32 bins histograms arc recorded. To calculate 
error bars they are transformed into jackknife bins along 
the lines of chapter 5.1 of^^. 

All calculations can be done by running one 2 GHz 
PC for about four weeks. As the runs at different pa- 
rameter values are independent, the real time is con- 
siderably shorter when several PCs are available. The 
Wang-Landau recursion^^ consumed never more then 
a few percent of a run. Cycling^^ and a flatness of 
-ffminZ-f^max > 0.5 was Considered sufRcient for iterating 
the Wang-Landau refinement factor. Such a crude flat- 
ness is sufficient when one docs not intend to converge 
into a reliable estimate of the spectral density, as origi- 
nally proposed in^^, but aims only at obtaining a working 
estimate of the MUCA weights. To use the Wang-Landau 
algorithm in this ways as a recursion for the first part of 
a MUCA simulation was suggested in Ref.^^. 

For the ho = 0.5 run on our largest lattice the MUCA 
energy histogram of the production part is shown in 
Fig. 3. The value Hq = 0.5 converts for this lattice to 
h = 0.500173 so that (Q)o = 1992 = K. In Fig. 4 
the energy histogram is restricted to entries for which 
Q = K = 1992 holds. It is this histogram, which is 
reweighted to P = and then normalized, so that its 
sum over energies, the denominator of the right-hand 
side of (9), matches the binomial coefficient (7). To 
monitor the entire (Q, E) distribution a histogram array 
^MUCA^g^ £;) of size N"^ would be needed. In our simula- 
tions we avoided arrays of size by focusing reweighting 
on one selected Q value, such that an array of size N is 
sufficient. However, this restriction to a microcanonical 
state was possibly not a wise decision. Compared to the 
analysis of^^ we find spurious fluctuations and increased 
error bars. Likely that could be smoothed out when the 
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FIG. 5. Overlap histogram from groundstates sampled. 



full array is available, which would for our largest lattice 

still fit into the memory of a PC, and in the analysis al- 
low to sum over Q for the normalization. As this would 
require to repeat all simulations, we cannot pursue this 

issue further at this point. 

The overlap histogram as measured in the correspond- 



ing groundstate distribution {Eg 



-2N 



-3200 for 



this lattice) is depicted in Fig. 5. Properly normalized 
the number of configuration per molecule follows from 
this histogram by using Eq. (9) for / values which are 
sampled with suSicient statistics. The cut-off values for 
sufficient statistics for / were determined from one half of 
the maximum value of H^^^^ = maxf[H^^'^^{f, Eg)] 
in the following way: 



/l < / < /2 



with 



A = mm [f;H^^^^{f,Eg) > , 



(11) 



(12) 
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FIG. 6. Wi(/) in the approximation (1) versus MUCA. 
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FIG. 8. (Color online) Difference between MUCA esti- 
mates and the approximation (1). The lines are only drawn 
to guide the eyes. 




FIG. 7. Wi{f ) in the approximation (1) versus MUCA. 



h = max [/; i/^^^^/, E,) > ^^.^^2] • (13) 



IV. ENTROPY ESTIMATES 

Fig. 6 compares the approximation Wi{f) of Eq. (1) 
with the estimates from our smallest lattice and Fig. 7 
with the estimates from our largest lattice. The differ- 
ences between the numerical results and the analytical 
approximation are in both cases small, but well outside 
the range of the numerical error bars. The latter point 
is demonstrated in Fig. 8, where we plot 



AW^if) = M^f N) - <(/) 



(14) 



for iV = 128 and 1600. A feature of Figs. 7 and 8 is that 
only patches of / are covered by the N — 1600 data. 
Each ho value defines such a patch by means of Eq. (11). 



TABLE II. Infinite volume extrapolations of Wi{f) and 
AWi{f) (the error bars of both quantities are the same). 



/ 


Wiif) 


AWiif) 


,/ 


Wiif) 


AWiif) 


0.50 


1.50620 (32) 


0.00620 


0.80 


1.27729 (26) 


0.00350 


0.65 


1.44166 (26) 


0.00587 


0.93 


1.09849 (20) 


0.00085 



The one corresponding to ho = 0.5 can be read off from 
Fig. 5: 0.7775 < / < 0.8075. By adding simulations for 
further ho values the uncovered / regions can be filled. 
We abstained from doing this, because it is only of aca- 
demic interest. Our corrections to the analytical approx- 
imation (1) show that this approximation is sufficiently 
accurate for practical applications, because error bars of 
experimental entropy estimates (e.g.,^*) are much larger 
than the correction to (1). 

Fig. 8 shows also the finite size corrections to 
W^^'^^{f,N) encountered when moving from N = 128 
to = 1600 molecules. These estimates together with 
those from the N = 360, 576 and 896 lattices allow 
one to perform infinite volume extrapolations Wi (/) = 
limAT^oo Wf^^^^if, N). As in Ref.^^ for the case / = 0.5 
we fit to the form 



l¥f^^^^(/,7V) = Wiif) + aN' 



(15) 



With the present data the 3-parameter fits turn out to be 
unstable and we reduce them to stable 2-parameter fits 
by using 9 = 0.92 from^"^ on input. For four / values the 
thus obtained infinite volume extrapolations Wi (/) are 
collected in table II (error bars are given in parenthesis). 

For / — 0.5 the fit is shown in Fig. 9. For the other 
/ values the shapes of the fits are quite similar. For 
/ — 0.93 the N — 128 estimate cannot be included, be- 
cause it would spoil the consistency of the fit. Remark- 
able is that finite size corrections for our microcanonically 
normalized data in Fig. 9 are much smaller than those in 
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FIG. 9. Finite size fit for Wi{N). 

the corresponding figure of Ref.^^, where a canonical nor- 
malization (summed over all Q values) of the density of 
states was used. Further the sign of the correction is 
opposite to that Ref.^^. As before, the present estimate 
is in good agreement with Nagle^*^, midershooting now 
his value slightly, whereas the value of^'^ is overshooting 
Nagle's estimate somewhat. 

The Wi (/) estimates together with their error bars are 
also plotted in Fig. 8. Besides for / = 0.93 they are only 
visible in the color version of this figure, because they 
fall within the error bars of the N = 1600 data. Inter- 
estingly the / = 0.93 extrapolation is considerably larger 
than the N — 1600 estimate and the sign of the correc- 
tion with respect to the approximation (1) flipped. 
While on all our lattices we have for sufficiently large / 
a crossover of the correction from positive to negative, 
this feature may disappear in the N ^ oo limit, so that 
the corrections are ultimately all positive. To illustrate 
lattice artifacts in the / ^ 1 limit we plot in Fig. 10 the 
AWi{f) values for / > 0.8. It is clear that the closest 
values to / = 1 reflect lattice artifacts and should not be 
used for the N ^ oo approximation. Still estimates for 
all values of / can be obtained, because the / range of 
the artifacts shrinks 1/N. 

V. SUMMARY AND CONCLUSIONS 

Our main finding is that the corrections to the analyti- 
cal approximation (1) are small. As illustrated in table II, 
they are never larger than Nagle's^" already small correc- 
tion to Pauling's^ value W^^*^"''"^ = 1.5. For the entropy 
this translates into 

AS < In (<^^' j - In (<-'-s) « « 0.46% . 

(16) 



< 




0.8 0.85 0.9 0.95 1 

f 



FIG. 10. Enlargement of the / -> 1 region for AWi. The 
lines are only drawn to guide the eyes. 

This is beyond the accuracy of nowadays measurements. 
But who knows about twenty years ahead? The veri- 
fication of the correctness of predicted correlations be- 
yond the Pauling-like approximation would be an ulti- 
mate confirmation of our understanding of ice. 

It is straightforward to include additional parameters 
in our approach, as introduced by the equations of Mac- 
Dowell et al.^*^. Each choice of parameters requires a 
simulational effort similar to that of Ref.^'^. So it would 
be tedious to map out corrections for the entire parame- 
ter space. In particular, we did not pursue this further, 
as due to our present results one may conjecture that 
these corrections are also small. If one likes to perform 
a check for a special choice of parameters, for instance 
because ongoing experimental measurements, the details 
given in our paper should allow researchers to set up the 
necessary simulations. 

Finally, there may well be applications of our approach 
to systems for which corrections to existing approxima- 
tions are not be small. For example, the method allows 
one to calculate the combinatorial entropy of small clus- 
ters of hydrogen bonds directly. They are observed as 
formation of ice layers in nanotubes^^ and expected to 
be of importance in the interaction of water with pep- 
tides, proteins and other biomolecules. Through a better 
understanding of their entropy insights derived from the 
study of ice may well lead to a better understanding of 
models, which have primarily been constructed to reflect 
interactions of water at room temperature (see'^^ for an 
overview). 
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